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Dynamical heterogeneities - strong fluctuations near the glass transition - are believed to be cru- 
cial to explain much of the glass transition phenomenology. One possible hypothesis for their origin 
is that they emerge from soft (Goldstone) modes associated with a broken continuous symmetry 
under time reparametrizations. To test this hypothesis, we use numerical simulation data from four 
glass-forming models to construct coarse grained observables that probe the dynamical heterogene- 
ity, and decompose the fluctuations of these observables into two transverse components associated 
with the postulated time-fluctuation soft modes and a longitudinal component unrelated to them. 
We find that as temperature is lowered and timescales are increased, the time reparametrization fluc- 
tuations become increasingly dominant, and that their correlation volumes grow together with the 
correlation volumes of the dynamical heterogeneities, while the correlation volumes for longitudinal 
fluctuations remain small. 
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For systems in the vicinity of the glass transition, ex- 
periments and simulations have shown the emergence of 
spatially heterogeneous dynamics (SHD): mesoscopic re- 
gions relax either much faster or much slower than neigh- 
boring regions 1-8 . SHD is believed to be crucial to the 
understanding of non-exponential relaxation, the break- 
down of the coupling between translational diffusion and 
viscosity, and even possibly the slowdown of the dynam- 
ics itself [U [5] . The origin of SHD is still uncertain, in 
part because of the lack of direct microscopic tests to at- 
tempt to disprove proposed theories [HJ IMTT) . Here we 
apply one such test [12] for the hypothesis that SHD is as- 
sociated with fluctuations in the time variable [TTHT5HT4] . 
and find that our molecular dynamics data are consis- 
tent with the hypothesis. This test can also be applied 
to particle tracking experimental data in colloidal 4 and 
granular systems [5], thus allowing to investigate a pos- 
sible unified explanation of SHD in diverse systems. Our 
results highlight that non-trivial correlation functions in 
the time domain contain useful information for the un- 
derstanding of SHD. 

As a glass-forming liquid approaches the glass transi- 
tion, its relaxation time and viscosity grow by many or- 
ders of magnitude, until the system can no longer equi- 
librate in laboratory timescales, i.e. it has entered the 
glass state pp. In equilibrium, the correlation function 
C(t,t w ) between the states of the system at the wait- 
ing time t w and the final time t depends only on t — t w , 
but if the system is out of equilibrium, it may display 
aging, i.e. a nontrivial dependence on both t and t w . 
Dynamical heterogeneity can be probed by defining a 



coarse grained local two-time correlation C r (£, t w ), which 
probes how much each individual region of the sample has 
changed between time t w and time t. "Fast regions" have 
small values of C T (t,t w ) and "slow regions" have values 
of C T (t,t w ) closer to 1. Thus the fluctuations of C r (t,t w ) 
represent the dynamical heterogeneity, and theories at- 
tempting to explain SHD should be able to explain those 
fluctuations. One of the proposed mechanisms for the 
origin of dynamical heterogeneity postulates that they 
are associated with local fluctuations in the time vari- 
able dUUMa], t ->• h r (t), i.e. 



C r {t,t w )=C(h T (t),h r (t w )), 



(1) 



where C(t,t w ) = C s } t, a \(t,t w ) is the global two-time cor- 
relation. This proposal originated in analytic calcula- 
tions in spin glass models in the long time limit that 
showed the presence of a broken continuous symmetry 
under reparametrizations of the time t — > hit) [i~3l I16j. 
which should give rise to the presence of Goldstone modes 
as described by Eq. [II Indirect evidence in favor of the 
presence of this kind of fluctuation in atomistic mod- 
els of glasses has been presented in [T51 [T7HTU] . In the 
present work, we introduce a more direct test, based on 
decomposing fluctuations into a transverse part satisfy- 
ing Eq. (fTl) and a longitudinal part containing all other 
fluctuations [5D] . This procedure allows one to separately 
quantify the strength and spatial correlations of both 
kinds of fluctuations, as a function of temperature and 
timescales, for a variety of glass-forming models, and is 
easily applicable to experimental data in glassy colloidal 
and granular systems. 



To probe fluctuations in structural glasses, we use [TS] 
C r (t,t w ) = iv(^) Er J (t l „)ei3 r C0S (q ' MO - r j(0))- 
Here ij(£) is the position of particle j at time t, B r de- 
notes a small coarse graining box around the point r, 
and the sum runs over the N(B r ) particles present in 
the coarse graining box at the waiting time t w . The 
global correlation function C(t,t w ), defined by extend- 
ing the average to all of the N particles in the system, 
is the self part of the intermediate scattering function. 
We have chosen the wavevector q to be at the main 
peak of the static structure factor S(q) for each sys- 
tem. We performed classical Molecular Dynamics sim- 
ulations of systems of N particles (1000 < N < 8000) 
that were equilibrated at high temperature Ti^> T g , then 
instantaneously quenched to a final temperature T and 
allowed to evolve for times several orders of magnitude 
longer than their typical vibrational times JT51 IT71 - H9| . 
We generated eight datasets by simulating four atom- 
istic glass-forming models [21]: an 80:20 mixture of parti- 
cles interacting via Lennard- Jones (LJ) potentials [T5ll2"2"] 
(dataset C), an 80:20 mixture of particles interacting 
via purely repulsive Weeks- Chandler- Andersen (WCA) 
potentials [T^j (datasets D-H), and short (10- monomer) 
polymer systems |17j interacting via either LJ potentials 
(dataset A) or via WCA potentials (dataset B). Near- 
est neighbors along the polymer chains are held together 
by FENE anharmonic spring potentials [17]. The ratio 
of the final temperature T to the Mode Coupling crit- 
ical temperature T c [23] was T/T c ~ 0.9 for datasets 
A-D, T/T c = 1.10 for datasets E-F and T/T c = 1.52 for 
datasets G-H. For datasets F and H, the samples were in 
equilibrium, but for all the others the samples were aging. 
Each dataset includes between 100 and 9000 independent 
runs with the same parameters. 

To test the hypothesis given by Eq. (HI) , we will use the 
fact that for our data l24l: 



now define $> a b,> 



=- I' 1 lCr(t a ,t b )], with a, 6 G {1,2,3}, 

$23,r/V $ 13,r and F r (*l, *2, *3) = 



C{t,t w )^f{h(t)/h(t w )). 



(2) 



where / can be fitted with a form such that C(t,t w ) re- 
duces to a stretched exponential in the equilibrium case: 
f(x) = q E Aexp{~[ln(x)/e } 13 } [25]. Here q E A, P and 9 
are fitting parameters that vary little from one dataset 
to another. However, the dependence of the a relax- 
ation time r on t w is quite different in the different sys- 
tems we consider [2B], and this leads to different forms 
for h(t) [2"3] : for aging polymers h(t) = exp[ln Q (t/i )], 
for aging particles h(t) = exp[(t/to) Q ], and in equilib- 
rium h(t) = exp(t/t ). We define $ Qb = f^ 1 [C(t a ,t b )}, 
with a, b e {1,2,3}. If Eq. ([2| is satisfied, we have 
&ab ~ h(t a )/h(tb), and we therefore obtain a triangu- 
lar relation [S] $13 w $i2 < i ) 23- In terms of the variables 
X = <&23/V$i3 and Y — 3>i2/V*i3j this l ea ds to the 
prediction that 1 ra XY, which is satisfied to a good 
approximation for all times and all of our datasets [25] . 
By using Eq. (f2j) , we now re-express our hypothesis, 
Eq. 0, in the form C T (t,t w ) = f (h r (t)/h r (t w )). We 



$i2,r/-v/^i3, r) whose fluctuations also encode the prop- 
erties of the dynamical heterogeneities. If the hypothesis 
in Eq. (IT]) is satisfied, then 1 = X r Y r , i.e., the relation 
holds locally not just globally. Since time reparametriza- 
tion symmetry is a long time asymptotic effect associated 
with glassy behavior, we expect that as the temperature 
becomes lower, the timescales become longer, and the 
system becomes more glassy, the probability distribution 
p(X r ,Y r ) should become anisotropic, and extend mostly 
in the direction of the global curve 1 = XY and not away 
from it. In other words, if we decompose the fluctua- 
tions representing the dynamical heterogeneity into lon- 
gitudinal and transverse variables [5D] , the fluctuations of 
the longitudinal variable er r = -4= m($i2,r$23, r/$i3,r) = 
-4= \n(X T Y r ) should become weaker than the fluctuations 
of the transverse variables ni r = 4= ln($i2, r /$23,r) = 
^ \n(Y r /X r ) and n 2r = i In(i> 12 , r $ 2 3,r$i3,r)- 

In Fig. ([I]) we show our results for p(X r ,Y r ). Because 
we are trying to detect collective fluctuations, we coarse 
grain over moderately large regions, containing on aver- 
age 125 particles. For each dataset, we find three triads of 
times t\ > t 2 > £3 such that {X{tx,t%,t'i),Y{t\,t%,t5J) w 
(0.8, 1.25), (1.00, 1.00), and (1.25,0.80) respectively. For 
each dataset and time triad, we show three contours of 
constant probability density p(X r ,Y r ), respectively en- 
closing 25%, 50% and 75% of the total probability. For 
datasets A-D, with T/T c ~ 0.9, the contours indeed fol- 
low the curve 1 = XY. This is more noticeable for the 
25% contour, which encloses the most likely fluctuations, 
than for the 50% and 75% contours, which additionally 
include rarer events. For datasets E and F, correspond- 
ing to T/T c = 1.1, the contours are still anisotropic and 
oriented along the direction of the global curve, but less 
so than in A-D, while for G and H, corresponding to 
T/T c = 1.5 the fluctuations away from the global curve 
are the strongest. For the higher temperatures, we find 
that the contours obtained in the aging regime (F, H) are 
similar to the ones obtained in the equilibrium regime (E, 
G) at the same temperatures Q]5]. These results can be 
directly connected to the fact that, as the temperature is 
increased, the separation of timescales is less pronounced, 
the finite time corrections to the time reparametrization 
symmetry become larger, and the effect of local time vari- 
able fluctuations become weaker. 

We now turn to a more quantitative analysis of the 
connection between the transverse fluctuating variables 
ttid 7r 2r, the longitudinal fluctuating variables a r , and the 
dynamical heterogeneity. A more detailed version of this 
analysis will be presented elsewhere [26] . Here we report 
results for fixed C(ti,ts) — 0.23, but similar results are 
obtained for other values of C(ii,i3) [55]. In the top 
panel of Fig. |2]) we show the ratio between the variances 
of the local transverse and longitudinal fluctuations as 
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FIG. 1: (Color online) 2D contours of constant joint probability density p(X r ,Y r ) = p($23,r/($i3,r) , $i2,r/($i3.rj 
computed using coarse graining boxes containing 125 particles on average. Each set of three concentric contours is chosen so 
that they enclose 25%, 50% and 75% of the total probability. Each panel from A to H contains results from the corresponding 
dataset, for (X,Y) ss (0.80, 1.25), (1.00, 1.00) and (1.25,0.80), with the times chosen as late as possible within each dataset. 
The global values (X(ti,t2,ts),Y(ti,t2,t3)), for all times ti > t% > £3 in each dataset, are shown with red points. 



a function of Xi,C(ti,t 3 ) 13 H8| 127] . which quantifies the 
strength of the dynamical heterogeneities. Similarly, in 
the bottom panel of Fig. (pi) we plot the ratio between 
the correlation volumes |27j of transverse and longitudi- 
nal fluctuations as a function of Xi,C(t±,ts)- I 11 both cases, 
we find that there is an anisotropy in favor of the trans- 
verse fluctuations, which grows as the strength of the 
dynamical heterogeneity increases. In particular, both 
ratios grow as the temperature is decreased, and in the 
case of systems in the aging regime, both ratios grow as 
the system relaxes, since X4,c(ti,t 3 ) is a growing function 
of t w at fixed C{t x ,H) P]-' 

Our hypothesis is that the dynamical heterogeneity 
originates in the Goldstone modes associated to fluc- 
tuations in the time reparametrization, as described by 
Eq. . We thus expect that the correlation length of the 
dynamical heterogeneity should be similar to the corre- 
lation lengths of the transverse variables 7Ti and 7T2, and 
that the longitudinal variable a should be short-range 
correlated. In Fig. pi), we show that this is indeed the 
case: the normalized correlation volumes [57] V ni , K- 2 
for the transverse fluctuations are approximately propor- 
tional to those for the dynamical heterogeneities, Vc 13 , 
and in particular they grow as the temperature is re- 
duced or as aging systems relax. By contrast, the nor- 
malized correlation volume for longitudinal fluctuations 
V a is essentially constant for all systems, temperatures 



and time regimes, and approximately equal to unity, in- 
dicating that the spatial correlations of the variable a do 
not extend beyond the coarse graining region. 

In conclusion, we have applied a stringent microscopic 
test for the hypothesis that dynamical heterogeneity in 
structural glasses is associated with the presence of spa- 
tially correlated fluctuations in the time variables, and 
we have found that all our results are consistent with 
this hypothesis. We have used data from molecular dy- 
namics simulations of atomistic systems to apply the test, 
but the same procedure can be applied to particle track- 
ing data from colloidal [4j and granular systems [8], and 
slight modifications would allow the study of light scat- 
tering [7] or dielectric noise [5] data. This opens the 
door to investigating the possibility of a unified theoret- 
ical explanation of dynamical heterogenity for molecu- 
lar liquids, colloidal liquids and granular systems. Our 
results highlight the advantages of studying dynamical 
heterogeneity by probing fluctuations of regions of the 
system, rather than probing individual particle fluctua- 
tions, since the latter will necessarily contain both collec- 
tive and non-collective components that are difficult to 
separate cleanly. They also highlight the fact that more 
complex correlations in the time domain contain informa- 
tion that is useful for the understanding of heterogeneous 
dynamical behavior. 
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FIG. 2: (Color online) Anisotropy ratios between trans- 
verse and longitudinal fluctuations, plotted as functions of 
X4,C(*i,t 3 ) = X4(<i,ts) [13 EH]! which measures the strength 
of the dynamical heterogeneity. Plotted for all datasets and 
all times consistent with C(ti,t^) — 0.23 and C(ti,t2) = 
C(t2,tz). Top panel: Ratio between the variances of the 
transverse and longitudinal fluctuations. Bottom panel: Ra- 
tio between the correlation volumes [37J of the transverse and 
longitudinal fluctuations. 
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FIG. 3: (Color online) Normalized correlation volumes V V1 
(crosses), Vn 2 (triangles), and V a (circles) for the transverse 
fluctuations 7Ti r , 7i~2r and the longitudinal fluctuations o T , plot- 
ted as functions of the normalized correlation volume Vc 13 = 
Vc(t ly t 3 ) of the dynamical heterogeneities [27] . Plotted for all 
datasets and all times consistent with C(ti,tz) — 0.23 and 
C{ti,ti) = C(t2,t3). The color key for the datasets is as in 
Fig-i 
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